単回帰 y=a+bx における,H1:a=0 と H2:b=0 の多重比較
data(thuesen, package = "ISwR")
d <- thuesen
str(d)
## 'data.frame': 24 obs. of 2 variables:
## $ blood.glucose : num 15.3 10.8 8.1 19.5 7.2 5.3 9.3 11.1 7.5 12.2 ...
## $ short.velocity: num 1.76 1.34 1.27 1.47 1.27 1.49 1.31 1.09 1.18 1.22 ...
res <- lm(short.velocity ~ blood.glucose, d)
ctab <- as.data.frame(summary(res)$coef)
round(ctab, 4)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.098 0.1175 9.345 0.0000
## blood.glucose 0.022 0.0105 2.101 0.0479
b <- coef(res)
V <- vcov(res)
m <- nrow(d) - length(b)
C <- diag(2)
s <- sqrt(diag(C %*% V %*% t(C)))
t <- C %*% b/s
R <- cov2cor(C %*% V %*% t(C))
# q : adjusted p value
library(mvtnorm)
q <- sapply(abs(t), function(x) 1 - pmvt(-c(x, x), c(x, x), corr = R, df = m))
ctab$adjust.p <- q
round(ctab, 4)
## Estimate Std. Error t value Pr(>|t|) adjust.p
## (Intercept) 1.098 0.1175 9.345 0.0000 0.0000
## blood.glucose 0.022 0.0105 2.101 0.0479 0.0637
# u : critical value for siginificant level 0.05
u <- qmvt(p = 0.95, tail = "both.tails", df = m, corr = R, delta = c(0, 0))$quantile
# CI for 95% family-wise confidence level
ctab$lwr <- C %*% b - u * s
ctab$upr <- C %*% b + u * s
round(ctab, 4)
## Estimate Std. Error t value Pr(>|t|) adjust.p lwr upr
## (Intercept) 1.098 0.1175 9.345 0.0000 0.0000 0.8367 1.3590
## blood.glucose 0.022 0.0105 2.101 0.0479 0.0637 -0.0013 0.0452